from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *
import SuperNeuron as SuperNeuronCl
import MyConstants as mc

lengthD = 637 * um
lengthS = 106 * um
lengthIS = 35 * um

radiusD = 2.34 * um
radiusS = 2.34 * um
radiusIS = 1.5 * um

areaD = lengthD * 2 * radiusD * pi
areaS = lengthS * 2 * radiusS * pi
areaIS = lengthIS * 2 * radiusIS * pi

ratioS = 1
ratioD = areaD/areaS
ratioIS = areaIS/areaS

Rin = 70 * ohm * cm

gleak = (1.0 / (3000.0 * ohm)) * cm ** -2

C = 1 * uF * cm ** -2

gnaS = 200.0 * msiemens * cm ** -2
gnaIS = 600.0 * msiemens * cm ** -2

gkS = 250 * msiemens * cm ** -2
gkIS = 200 * msiemens * cm ** -2

ENa = 115*mV
EK  = -5*mV


rcDtoS = (Rin * (lengthD / 2 + lengthS / 2)) / (radiusD ** 2 * pi)
gcDtoS = (1 / rcDtoS)/areaS

rcStoIS = (Rin * (lengthS / 2)) / (radiusS ** 2 * pi) + (Rin * (lengthIS / 2)) / (radiusIS ** 2 * pi)
gcStoIS = (1 / rcStoIS)/areaS

tausyn1 = 1*ms
tausyn2 = 20*ms

Ee = 60*mV
Ei = -15*mV

class RNeuronFactory (SuperNeuronCl.SuperNeuron):
    eqs=Equations('''alpha_mS = 1/ms*(0.8*(22.5-vm/mV)/(exp((22.5-vm/mV)/4)-1)) : 1/ms
        beta_mS = 1/ms*(0.7*(vm/mV-50)/(exp((vm/mV-50)/5)-1)) : 1/ms
        alpha_hS = 1/ms*(0.32*exp((35-vm/mV)/18)) : 1/ms
        beta_hS = 1/ms*(10/(exp((50-vm/mV)/5)+1)) : 1/ms
        alpha_nS = 1/ms*(0.03*(25-vm/mV)/(exp((25-vm/mV)/5)-1)) : 1/ms
        beta_nS = 1/ms*(0.5*exp((20-vm/mV)/40)) : 1/ms
        
        alpha_mIS = 1/ms*(0.8*(15-vIS/mV)/(exp((15-vIS/mV)/4)-1)) : 1/ms
        beta_mIS = 1/ms*(0.7*(vIS/mV-40)/(exp((vIS/mV-40)/5)-1)) : 1/ms
        alpha_hIS = 1/ms*(0.32*exp((40-vIS/mV)/18)) : 1/ms
        beta_hIS = 1/ms*(10/(exp((40-vIS/mV)/5)+1)) : 1/ms
        alpha_nIS = 1/ms*(0.03*(15-vIS/mV)/(exp((15-vIS/mV)/5)-1)) : 1/ms
        beta_nIS = 1/ms*(0.5*exp((10-vIS/mV)/40)) : 1/ms
        
        dmS/dt = alpha_mS*(1-mS)-beta_mS*mS : 1
        dhS/dt = alpha_hS*(1-hS)-beta_hS*hS : 1
        dnSo/dt = alpha_nS*(1-nSo)-beta_nS*nSo : 1
        
        dmIS/dt = alpha_mIS*(1-mIS)-beta_mIS*mIS : 1
        dhIS/dt = alpha_hIS*(1-hIS)-beta_hIS*hIS : 1
        dnIS/dt = alpha_nIS*(1-nIS)-beta_nIS*nIS : 1
        
        ISynI = ipsp*(vD-Ei) : amp*umetre**-2
        ISynE = epsp*(vD-Ee) : amp*umetre**-2
        
        dvm/dt = -(gnaS*mS**3*hS*(vm-ENa) + gkS*nSo**4*(vm-EK) + gleak*(vm-ElS) + gcDtoS/ratioS*(vm-vD) + gcStoIS/ratioS*(vm-vIS))/C : mV
        dvIS/dt = -(gnaIS*mIS**3*hIS*(vIS-ENa) + gkIS*nIS**4*(vIS-EK) + gleak*(vIS-ElD) + gcStoIS/ratioIS*(vIS-vm))/C : mV
        dvD/dt = -(ISynI+ISynE+gleak*vD+gcDtoS/ratioD*(vD-vm))/C : mV
        
        ElS : mV
        ElD : mV
        ''')
        
        
    def init(self, size,synapse='exp',Evariable=1,rand=False):
        self.size = size
        self.rand = rand
        self.Evariable = Evariable
        self.add_synapses(synapse,tausyn1,tausyn2)
        self.add_rand_curr(rand)
        
        self.Group=NeuronGroup(self.size,model=self.eqs,
                      threshold=EmpiricalThreshold(threshold=40*mV,refractory=5*ms),
                      implicit=True, clock=Clock(dt=0.01 * ms))        
        self.reinit()
       
        
        return self.Group  
    def reinit(self):
        self.Group.vm = 0 * mV
        self.Group.vD = 0 * mV
        self.Group.vIS = 0 * mV
        self.Group.mS = 0
        self.Group.hS = 1
        self.Group.nSo = 0
        self.Group.mIS = 0
        self.Group.hIS = 1
        self.Group.nIS = 0
        self.Group.ElS = randn(self.size)*2*self.Evariable*mV
        self.Group.ElD = randn(self.size)*2*self.Evariable*mV
        
        